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A new finite volume (FV) discretisation method for the Lattice Boltzmann (LB) equation which 
combines high accuracy with limited computational cost is presented. In order to assess the per¬ 
formance of the FV method we carry out a systematic comparison, focused on accuracy and com¬ 
putational performances, with the standard streaming (ST) Lattice Boltzmann equation algorithm. 

To our knowledge such a systematic comparison has never been previously reported. In particular 
we aim at clarifying whether and in which conditions the proposed algorithm, and more generally 
any FV algorithm, can be taken as the method of choice in fluid-dynamics LB simulations. For 
this reason the comparative analysis is further extended to the case of realistic flows, in particular 
thermally driven flows in turbulent conditions. We report the first successful simulation of high- 
Rayleigh number convective flow performed by a Lattice Boltzmann FV based algorithm with wall 
grid refinement. 

Keywords: Lattice Boltzmann Method, Finite Volume approach, Rayleigh-Benard system, turbulent convec¬ 
tion 


I. INTRODUCTION 

Since its introduction, the Lattice Boltzmann (LB) equa¬ 
tion method for fluid-dynamics simulations has enjoyed 
increasing success !l]]. Reasons are two-fold, on one 
hand the mesoscopic level of description on which it is 
based goes beyond the Navier-Stokes continuum mat¬ 
ter description of fluids and it eases, as compared to 
other macroscopic methods, to accommodate for com¬ 
plex effects such as for instance the interaction between 
different fluid components, phase-change processes, non- 
newtonian rheology. The extensions of the LB meth¬ 
ods in such directions are countless ( e.g. multiphase and 
multicomponent flows 0 - 0 ; flows with suspensions 0 
fl2| ; emulsions [l3| ; porous media JTM3 ; natural convec¬ 
tion (HI; reactive transport fsl |20H : combustion [U [22j ; 
magneto-hydrodynamics [23|, [2J]). On the other hand, 
the LB method has also very appealing features from a 
computational point of view. It is simple to implement, 
free of numerical diffusion and stability issues and suit¬ 
able for parallelization due to its local-in-space character. 
However, when it comes to the simulation of turbulent 
flows one shortcoming of the method, the limitation to 
equispaced grids , becomes evident. We recall here that a 
developed turbulent flow in the presence of any sort of 
bounding geometry (or any local forcing term) develops 
space inhomogeneities and as such grid refinement in nu¬ 
merics becomes necessary. It shall be made clear that 
in such a context, grid-refinement is not an additional 
requirement in order to increase the accuracy of a simu¬ 
lation but is rather an unavoidable need in order to save 
memory usage and computational power and being able 
to access higher - read more realistic - turbulent flows 
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regimes. In summary any state-of-the-art computational 
fluid dynamics method (CFD) calls for grid refinement. 
Several approaches have been proposed in order to over¬ 
come the shortcoming of equi-spaced grids in Lattice 
Boltzmann (LB) equation simulations. Here we mention, 
i) the grid refinement methods which make use of locally 
nested equispaced grids 05|, ii) the techniques based on 
off-lattice interpolation schemes [26l [27| iii) the finite- 
difference [28j . finite-volume §3 or finite-elements 

10 LB discretisation methods, and iv) the extension of 
the LB equation to general manifolds (32j. There are 
however important drawbacks, all such reformulation are 
computationally more expensive, or introduce extra sta¬ 
bility limitations enforced by space/time discretisation 
which were not present in the original so-called stream¬ 
ing based implementation. Presently, the only viable 
way seems to be the nested monospaced grid method (i), 
which has allowed to simulate turbulent channel flows 
[HI and even more complex flow geometries 0]. How¬ 
ever, in this case the advantage in terms of accuracy and 
efficiency compared to state of the art direct numerical 
simulation (DNS) e.g. spectral methods is limited. 

This paper presents a new Finite-Volume (FV) discreti¬ 
sation method for the Lattice Boltzmann equation which 
besides a high level of accuracy also displays a contained 
computational cost. In order to assess the performance 
of this new FV method we carry out a systematic com¬ 
parison with the standard streaming (ST) formulation. 
To our knowledge such a methodical comparison of accu¬ 
racy and computational performances has not been made 
previously. We aim at clarifying whether and in which 
conditions the proposed FV algorithm can be taken as 
the method of choice in fluid-dynamics LB simulations. 
The paper is organised as follow. In the next section we 
describe two different discretizations of the LB equation. 
To begin with, we briefly review the key points of the 
streaming implementation ( fil All , then we present the 
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new FV based formulation 1 HIIBI) . Our guidelines in the 
development of the new FV method are the the simplicity 
and computational efficiency of the implementation, yet 
retaining a level of accuracy which takes the ST method 
as the baseline. Results of this study are reported in 
m First, we address the accuracy comparison, later on 
the computational efficiency of the algorithm. A further 
section fi fVl) takes the analysis to more realistic flows, 
in particular thermal flows in turbulent conditions. To 
our knowledge, for the first time a high-Raylcigh number 
convective flow with wall grid refinement is simulated by 
a LB based code. Final remarks and perspectives are 
reported in the conclusions. 


pu = Ti a c a f a and the kinematic viscosity as v = r 
where the constant c s stands for the so-called lattice 
speed of sound, whose value depends on the specific ve¬ 
locity lattice topology (see for instance Succi’s book QJ 
for details). Finally, F a is a forcing term, constructed 
in such a way to model the effect of a macroscopic body 
force term. 

Note that eq. © is discrete just in the velocity space ( 
for this reason it is also known as discrete velocity Boltz¬ 
mann equation). Up to this level no discretisation has 
been taken neither in the spatial domain or in the tem¬ 
poral one. Such further discretisations can take different 
paths as we describe in the following sections. 


II. METHOD 

We focus here on the LB method with the Bhatnagar- 
Gross-Krook (BGK) collision operator, which is charac¬ 
terised by a single relaxation time r towards a local equi¬ 
librium state. The equation of motion reads: 

Bf 1 

-£+c a -Vf a = -(f?-f a )+F a a = 0,..., N pop (1) 

at t 

where f a (x, t) is one of the N pop distribution functions for 
particles (also called populations ) with velocity c a at po¬ 
sition x and time t. The set of f a distribution functions 
relax towards a local equilibrium state f^ q (x,t) which 
is prescribed in terms of local macroscopic variables (in 
non thermal models, as here, they are just the fluid ve¬ 
locity and density) [l|. The macroscopic fluid mass and 
momentum density can be computed as p = T, a f a and 

I 

At 

fa(x + C a At,t + At) = f a (x,t) + — 

T 


A. Outline of the Streaming Lattice Boltzmann 
algorithm 

It is here useful to briefly recall the steps that have to be 
made in order to obtain from eq. © the standard stream¬ 
ing (ST) LB algorithm. First, by applying to the above 
partial differential equation (PDE) the technique of the 
characteristics along the lines x(t) = cc(0) + c a t , one 
gets the ordinary differential equation (ODE) 

-U = -(f e a q -U) + F a ■ ( 2 ) 

at t 

Second, the discrete integration in time, of step At, is 
performed by applying the semi-implicit Crank-Nicolson 
method. Such a step is followed by a convenient redefi¬ 
nition of the distribution functions for the lattice popu¬ 
lations f a = fa — ^i(fa 9 — .fa + T~F a ) which makes the 
scheme explicit, leading to [351 ] : 


~^{ x ,t)-f a {x,t)) + At(l-^jF a (3) 


where r = r + At/2 is a redefined relaxation time 
(f > At/2). It is easy to derive the relations between 
the macroscopic variables and the tilded (~) quantities, 
they are respectively: p = T, a f a , pu = Y, a c a (f a + 4 fF a ) 
and is = (f — At/2)Cg. Note that a factor 1 — in front 
of the forcing term needs to be introduced a posteriori in 
order for the discretised eq.© to give the same hydro¬ 
dynamics limits as © M- 

The numerical implementation of © is straightforward. 
It can be divided in two steps: i) the computation of the 
right-hand-side and ii) the displacement (or streaming) 
of the computed values on the lattice according to the di¬ 
rection and intensity of c a . It is important to note that 
the integration along the characteristics introduces a link 
between the space and time discretization, which reads 
c a = Ax a /At. This means that if one choose the carte¬ 
sian components of the set of c a velocities to be either 
±1 or 0, it implies that Ax a ,i = At or 0. The standard 


choice (but not the only possible one) is Ax a i = At = 1 

0 . 


B. Lattice Boltzmann Finite Volume Formulation 

The method of characteristics is very convenient from a 
computational point of view because it reduces the com¬ 
plexity of the integration of a PDE to a simple ODE, 
however at the same time it introduces a tight link be¬ 
tween the shape of the velocity lattice and the spatial dis¬ 
cretisation mesh. Such a constraint can be removed if 
one takes the more usual numerical approach based on i) 
a direct spatial discretisation of eq. © combined with ii) 
an independent time discretisation phase. For the first 
step, several standard options are available, such as finite 
elements, finite differences or finite volumes methods. 
The idea of using a finite volume method to decouple the 
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spatial numerical mesh from the velocity lattice struc¬ 
ture was first proposed by Nannelli and S. Sued, J30] 
(see also |2l|), in this seminal paper a low-order upwind 
scheme was suggested for the discretisation of the ad- 
vection (or flux) term. The idea was further refined in 
Amati et al. 1371 ] . where piece-wise linear interpolation 
scheme was suggested for the treatment of the flux term. 
While these first works were limited to stretched carte¬ 
sian grids, Chen [38] presented a volumetric formulation, 
based on a cell-centered discretisation scheme, which al¬ 
lowed for the adoption of arbitrary structured meshes. 
The formulation was further developed by Peng et al. 
[HI, |4(j through cell-vertex FV scheme, which displayed 
enhanced stability properties. Sbragaglia & Sugiyama 
Ej applied Peng’s scheme to an energy-conserving LB 
model to study for the first time thermal convective flows. 
More recently Ubertini et al. [42| addressed the problem 
of unstructured bidimensional triangular meshes, which 
allow great flexibility on one hand, but also reintroduce 
known issues related to numerical stability. This has been 
further refined in a work by Zarghami et al. [43j through 
a cell-centered FV approach on arbitrary mesh in two di¬ 
mension. A Total Variation Diminishing (TVD) formu¬ 
lations for LB FV algorithm has been suggested by Patil 
et al. [HI], where stability and accuracy can be efficiently 
enhanced. Despite all these contributions, at present the 
situation is still far from being solved. If on one side it 
has been shown that a satisfactory level of precision can 
be reached by the FV method on the other hand this is 
often at the price of the high computational costs needed 


to obtain a stable algorithm. As an example in a recent 
work [43l l, where a series of laminar but relatively com¬ 
plex flows over non-homogeneous meshes were simulated, 
a fifth-order Runge-Kutta scheme had to be adopted for 
the time discretisation in order to have stable results. 
The consequence on the computational cost is evident 
since in such a scheme the advection terms of 0 need 
to be computed five times per time step (while in the ST 
method it is performed just by means of a memory shift, 
the streaming). As a consequence the FV LB is rarely a 
method of choice in fluid-dynamics simulations (see also 
the discussions in [45l] and in HD- 

The present paper further develops the finite volume 
Lattice Boltzmann method in order to simulate fluid flow 
problems with higher accuracy, greater stability prop¬ 
erties and comparable performance as the ST method. 
The FV method that we propose is of the type denoted 
as cell-centered (as opposed to vertex-centered, see Fig. 
COE- Its most original features concern the approach 
taken for the time discretisation fl jll B 21) and the method 
of fluxes computation which adopt, for the first time 
in this context, a quadratic upwind scheme f HII B 31) . 
In the following we detail the steps taken in developing it. 

1. Space discretisation 


Upon integration of eq. 0 over a volume V (of surface 
S) and by applying the flux theorem we get 


f v d ~^ dV + J s Ca ' nfa dS 


We then assume that every term in the volume integrals 
can be considered as constant and its magnitude taken 
at a reference location x (also called node) inside V. 

The term in the surface integral however, carries some 
kind of spatial variability. When such a surface is de¬ 
composed in M faces (as in a structured grid of nodes 
with connectivity index M) it is convenient to make the 
assumption that f a is constant on each of the Sj sur¬ 
faces perpendicular to rij and denoting its value with 


f -(/a 9 - /a) dV + [ F a dV (4) 

Jv T Jv 


[/a] .. This altogether leads to: 

+ f c a • nj [f a ] • = i(/^ - fot) + F a (5) 

Where summation over the repeated index j is ap¬ 
plied. 

2. Time discretisation 


If the time derivative is discretized by the explicit Euler 
scheme, we get: 


/, 


(i+At) _ At) 


.Si 


= - A tyC a 


Tin 


+-(f e a qW -f^) + AtF a , 


At 


reqr(t) 


(*)'! 


( 6 ) 


where the superscript indexes ( t ) and (t + At) denotes 
respectively the current and the next discrete time in¬ 


stant. Such an approach however, puts tight bounds on 
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the maximum allowed At. It is easy to show that if we 
discard the advection and the forcing terms and assume 
f eq to be constant, the stability region of the method is 
0 < At < 2 t. Empirically it is possible to show that this 
range becomes even narrower when the non-local advec¬ 
tion term, the forcing and the time dependency in f eq are 
taken into account. The fact that A t max depends on and 
is bounded by the value of r is a known problem in FV 
LB implementations. It poses, among others, a severe 
limitation for the simulations of turbulent flows ( i.e . low 
viscosity flows). On the opposite, such a constraint does 
not exist in the ST approach (where At is independent of 
r). Different solutions have been proposed in the litera¬ 
ture, often resorting explicit time discretisation schemes 
of higher order, for example multi-stages Runge-Kutta 
schemes. However, as we mentioned above such schemes 


only produce marginal improvements at the expenses of 
considerably increasing the computations. The Runge- 
Kutta schemes for example requires multiple evaluations 
of the full right-hand-side terms on (0. We opt for a 
different approach, with a better trade-off between the 
enhancement of the stability limit for At and the growth 
in computational cost. 

Similarly to what is done for the classic LB streaming 
method, in the steps from eq.0 to 0, a possible im¬ 
provement consists in taking also for the FV algorithm 
a semi-implicit integration scheme. However, this is not 
directly possible for eq. 0 because of the presence of the 
advection term. Therefore, we propose to limit such an 
approach only to the collision and forcing terms, after 
few manipulations (more details in the appendix) one 
gets the discretised form: 


fi t+At) = fa - At 


V 


fa + - fa) + 


+ - fa) + At fl ~^\F n 


(7) 


The above equation share the same definitions of 0 for 
the tilded distribution function, f a and the relaxation 
time (f). The rule of computing the macroscopic fields 
and the viscosity v = t <? s = (f — At/2) are exactly the 
same as for the ST algorithm. Correspondingly, the term 
1 — in front of the forcing has been introduced a pos¬ 
teriori to keep the same hydrodynamic limit. However, 
one can immediately note the advected field in the equa¬ 
tion is not simply a distribution function but a rather a 
complex term involving also the equilibrium distribution 
and the forcing. The main advantage of this approach 
is that a stability analysis under the same hypothesis 
mentioned above (neglecting advection, forcing and time 


dependences in the equilibrium function) shows now that 
every time step length At is stable. However, we find that 
the situation reached so far is not yet satisfactory. From 
simple numerical tests we observe that even with this 
discretisation scheme the time-step size is still restricted 
by the relaxation time, particularly for small relaxation 
time values. The origin of this still limited stability of the 
scheme lies now in the advection term. For this reason 
a further refinement is proposed. We set it into place by 
applying the so-called Heun predictor-corrector scheme 
to the advection term. In other words we use the calcu¬ 
lation of the population based on 0, now called /*, as 
an intermediate value for constructing an explicit trape¬ 
zoidal integration rule applied to the advection: 


f(t + A t)=f ~ a _ At ^ Ca 


f* _j_ ZAt / feq* f*\ I At rp 

J a ' 9 r & Jot/' 2 r ( 


fa + mfa q - fa) + f F t 


Tin 


J J 




( 8 ) 


where F* indicates the LB forcing term computed from 
/*. This scheme enjoys greater stability at the addi¬ 
tional computational price of a second evaluation of the 
advection term. In order to make this observations more 
quantitative we should first specify the way in which the 
flux terms [.. .]j are computed. Indeed, the exact sta¬ 
bility properties of the method depends upon the imple¬ 
mentation of the advection term, that we discuss in the 
following section. 


3. Approximation of the advection term 

There exist several ways to estimate the non-local term 
[f a ]j and each one can be characterised by a spatial order 
of accuracy. The complexity of such an estimation also 
depends on the grid characteristics. Even for structured 
but irregular grids an high-order estimation of [/„]. be¬ 
comes expensive in computation terms. In order to sim¬ 
plify such a problem, we limit the following discussion 
to the case of structured regular grids, that is to say to 
the case where the nodes lie on lines. This is the case 
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FIG. 1. (a) Cartoon of the finite volume space discretisation: 
the dot denotes the position of the cell center (where the value 
of f a {x) is defined), while the lines marks the cell boundaries. 
The cell has volume V and each boundary surface is denoted 
with Sj with j — 0,.... 3 in two dimensional space, (b) Sketch 
of the quadratic upwind interpolation scheme (QUICK) for 
the estimating the value of f a at the cell boundary position 
xs■ Note that the interpolation method make use of different 
nodes according to the direction of the population velocity 

Ca. 


for instance of a non-uniform cartesian grid (the typical 


case of wall-refinement), but it also apply to a uniformly 
skewed non-orthogonal grids. 

It has been long known that fluxes in advection equa¬ 
tions are better approximated by upwind schemes, which 
are interpolation schemes biased in the direction deter¬ 
mined by the sign of the characteristic speeds (the set 
of c a in our case). At the lowest order of accuracy, and 
easiest level of implementation, there exist the first order 
up-wind scheme, increasing the refinement leads to linear 
interpolation schemes or even to more refined quadratic 
schemes (which are of 3 rd order of spatial accuracy). 
While low-order schemes introduces artificial numerical 
dissipations, higher-order ones lead to spurious oscilla¬ 
tions, especially evident near the boundaries. This is 
also true in the present cell-centered FV implementation, 
in particular zero-order or linear up-wind interpolation 
schemes leads to inaccurate results. Even a cell-centered 
symmetric schemes, which here does not display extra 
dissipation, produces inaccurate results in the presence 
of boundaries. Empirically, we find that the quadratic 
upstream interpolation, known as QUICK method 1471] . 
is the simplest one to give accurate results both in open 
(be. periodic) and bounded domains. 

According to this approach, on each surface Sj at posi¬ 
tion say xs j , [f a \j is approximated via a combination of 
the value of f a in the two nodes bracketing the surface 
(denoted with x and x + ) and a third node that is located 
upstream respect to direction of the projection of c a on 
hj (denoted either x ++ or x~ ). The interpolant function 
is a parabola a + b £ + c £ 2 , with £ the linear coordinate 
spanning on the line connecting the nodes (see sketch in 
Fig. QJ>). This leads to: 


_I 

[faixs^j = (1 - 71 + 72 )f a (x) + 7l f a (x+) - 72 f a {x~) \ a . c a .fij >0 

+ (1 - 73 + 74 )fa(x + ) + 73 f a (X) - 74 fi(x ++ ) \ a , c {x -fij < Q (9) 


where the 71 , 2 , 3,4 are 4 coefficients, that shall be eval¬ 
uated and/or stored for each surface of the control vol¬ 
umes. 


4- Force term 

Finally a brief remark on the forcing term F a in the LB 
equation. The simplest way to implement it, is by the 
expression: 

r, _ c a ■ p a nn\ 

F a — „ (19) 


correct macroscopic effect of a body force term. However, 
when the body force is time/space dependent and eq. dTj) 
is discretised in space and time, such as in ([3j), the above 
expression needs to be refined in order to remove spu¬ 
rious discretisation terms that would otherwise appear 
in the macroscopic limit. The corrected expression, first 
proposed by Guo et al. jUJ, is 


F a — iCq. 




( 11 ) 


where the summation over index a is not implied, w a is a 
lattice dependent weight and the product pa represents 
the force per unit volume in physical space (for exam¬ 
ple in case of a gravitational external field a = g). The 
above expression satisfies the conditions T, a F a = 0 and 
Tj a c a F a = p a, which are required for eq. © to give the 


with an overall multiplicative factor 1 —At/(2r) when the 
distribution functions f a are evolved in place of f a . Note 
that accordingly (by employing the relation between f a 
and pu and p given in ] jH Al) one gets the fluid velocity 
as u = E„c a / a /E a / a + ^ a. 
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5. Boundary conditions 


In the following we consider the implementation two 
types of boundary conditions (BC): i) no-slip walls and 
ii) fixed density (or equivalently pressure) boundaries. 
The physical domain boundaries lie on the faces of the 
external control volumes. Similarly to the bounce-back 
approach for the streaming LB algorithm, we introduce 
in-wall ghost cells. However, in the QUICK treatment 
of the advection two ghost cells are needed instead of 
one. The ghosts cells are located in-wall and have cen¬ 
tres at position mirroring the first and second nodes in 
the fluid domain. Let’s suppose that the quantity to 
be advected is f a and that the boundary condition is 
to be imposed on the S cell surface, whose center is at 
xs . For simplicity we assume that S lies along the plane 
(■ y 7 z) perpendicular to x , with the x axis pointing in¬ 
ward ( i.e. in the fluid bulk direction). Consequently, the 
first two nodes in the fluid domain are located at position 
x\ = xs+Axi/2 and x 2 = xs+Ax\-\-Ax 2 /2 , where Axi 
and Ax 2 represents the linear size of the two first discreti¬ 
sation volumes. Accordingly, the ghosts cells are at posi¬ 
tions x_i = xg — Aaq/2 and x_ 2 = xs — Ax± — Ax 2 /2. 
A no-slip boundary condition requires u(xs ) = 0, while 
the density at p(x$) is free to take any arbitrary value. 
This corresponds to the constraint E a c a f a (x s ') = 0. The 
simplest way (but not the only one) to enforce it, is to 
set the in-wall nodes as the following: 

fa(X — l) finv(a)('E 1)7 

fa{x~ 2 ) /mu(a)(^ 2 )i (12) 

where inv(a) is an integer valued function that selects 
the population moving along the opposite direction with 
respect to c a . We have verified that such a choice does 
not introduce artificial fluctuations at the boundary, that 
would quickly generates instabilities. This implementa¬ 
tion of BC, that we dub double reflection , has a first-order 
of accuracy in space (it does not implement the quadratic 
interpolation) and therefore it leaves room for further im¬ 
provements. 

If instead we are interested to impose a density value at 
the border, say ps = E a fot{xs), we need to resort an 
extrapolation strategy. We proceed as follow, first the 
density value p(x_ i) is linearly extrapolated from the 
values ps and p(x i), similarly p(x_ 2 ) is derived from ps 
and p(x 2 ). Second, we assign the in-wall the distribution 
functions as follows 

/<*(*- 1 ) = f a (xi), 

fa(x_ 2 ) = f a (x 2 ). (13) 

P{x 2 ) 

Also the above choice, a rescaling of the bulk distribution 
functions, is not the only viable way for the implemen¬ 
tation of fixed density BC, however it is one that has 
revealed to not to introduce wall disturbances. As a final 
remark, we shall note that in the implementation of ([SI) 
the boundary conditions need not to be implemented on 
the redefined distribution function f a but rather on the 


original f a = f a + fy(/® 9 ~ fa) + ^r F a- 
The algorithm presentation given so far is independent 
of the particular microscopic velocity lattice topology. In 
the present work and for the accuracy study presented in 
the reminder of this manuscript we make the choice to al¬ 
ways use the so called D3Q19 lattice, which is a standard 
option for three-dimensional LB simulations and reduces 
to the D2Q9 lattice for two-dimensional flow problems 

PI- 


. „ born 
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FIG. 2. Illustration of the finite-volume arrangement for the 
implementation of the double-reflection boundary conditions. 


III. ACCURACY TESTS 

In this section we address the accuracy of the present LB 
FV algorithm and we compare it with the ST algorithm. 
In particular we approach the following questions: i) to 
which degree the FV algorithm correctly describes the 
dynamics of a low Reynolds number viscous flow? Which 
is its order of spatial accuracy and how does it compare 
with ST? ii) Is there any optimal usage of the FV algo¬ 
rithm in order to take advantage of the grid refinement 
and obtaining highly accurate solutions? 


A. Viscosity evaluation 

A simulation is performed on a physical domain of size, 
[L Xl L y , L z \ = [1, 64,1], For this test, the number of grid 
nodes per direction (indicated with N XtVtZ ) is also the 
same. The flow is initialised with an one-dimensional 
sinusoidal velocity amplitude profile of the form 

u(z, y, z) = ( u x (y), u y , u z ) = sin (^j^j , 0, 

(14) 

and it is left to decay in time. We monitor the be¬ 
haviour of the total kinetic energy in time, k to t(t), 
which is expected to decrease exponentially as k to t(t) = 
jA 2 Ly e - 2 ( 27r /W) v t , with v representing the fluid kine¬ 
matic viscosity. The reproduced value of v can be de¬ 
duced from a least-square fit of k to t{t) and then com¬ 
pared to the theoretically expected value v = r c 2 = 
(f — At/2) c 2 . The degree of accuracy of the FV method 
measured in such a way is compared with the ST method 
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and reported in figure [3] While it is known and expected 
that accuracy carries some form of dependency with the 
relaxation time r, we observe that both FV and ST meth¬ 
ods reach the maximal accuracy (minimal value of the 
error denoted E r ) around r = 0.5. Such behaviour has 
been already reported in the works of Holdych et al. (H} 
and Kruger et al. [49]], however ST in that very same 
case performs better by a factor 10. Moreover, in general 
the ST error grows less than the FV one for all r < 0.5. 



FIG. 3. Relative error of measured kinematic viscosity u num 
respect to the expected one Uth = r c? s as a function of the 
relaxation time r. In the finite volume case At = 1 for r > 
0.13 (marked with a vertical line) and At = 0.1 for r < 0.13, 
while in the Streaming case At = 1 always. In the inset, the 
absolute value of the same error in log-log scale. 


B. Steady Poiseuille flow 

Our second tests addresses the case of a simple bounded 
flow in the same spatial domain as above, [L x ,L y ,L z ] = 
[1, 64,1], The flow is initiated with a parabolic Poiseuille 
velocity profile U x (y) = 4 t/ max L“ 2 y (L y - y) corre¬ 
sponding to a Reynolds number Re = LyUmax/v = 10. 
A uniform volume forcing along the x direction and no- 
slip boundary conditions at y = 0 and y = L y positions 
are imposed, while periodicity is implemented along the 
horizontal direction, x. The forcing is implemented via a 
constant acceleration a x = -jp j U using equation dill) . 

The simulated flow profile (denoted with u x {y)) keeps the 
original theoretical shape U x (y) with tiny adjustments 
depending on the method. In order to compare these two 
functions we use the relative difference \\u x — XJ \\2 / ||£7||2 
where || ... ||2 denotes the L 2 norm, which in its discre- 


tised form is computed as: 


ft \ 1 ! 2 

\\f{x)h = yjjixfdxj = (T,* =1 f?Axi) 1/2 (15) 



FIG. 4. Relative error on the Re = 10 velocity Poiseuille 
flow profile at changing the number of grid points ( N ) and 
keeping fixed the grid spacing A* = L y /N y = 1. Proof that 
FV has same order of incompressibility accuracy of ST but is 
less accurate of a factor 8 to 10. 


In figure HI we show the L 2 relative difference results 
at varying the physical domain size in y direction, i.e., 
changing L y and at the same time N y (or in other words 
keeping fixed the grid spacing Ax = L y /N y = 1). Note 
that in this test we are actually varying the maximal 
Mach, Ma = U m a X /c s , number of the flow. Since the 
LB equation is 0(Ma 2 ) accurate Q, the figure proves 
that both the FV and ST methods posses the same level 
of incompressibility accuracy. However, we can clearly 
notice that ST in this conditions is still on average more 
accurate by a factor ~ 8 — 10 as compared to FV, such 
a difference is inherited from the behaviour at r = 0.5 
of the viscosity accuracy highlighted in the previous test. 

As further step, we address the effect of a stretched spa¬ 
tial grid on the overall accuracy of the Poiseuille flow 
simulation. To this end we implement three types of 
commonly used wall-normal stretched grids. The y- 
coordinate value of the cell volume centres (or simply 
nodes) is given by 

Vi = ^ +1 2 + ^ with 0 < N y 

where the the coordinates at the volume boundaries, 
are defined as 
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Chebychev nodes:f;= 
hyperbolic tangent :£*= 

hyperbolic sine:£i= 
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L and IV denote here respectively the physical domain 
size and the grid size (the sub-script index y have been 
dropped for brevity), and si, s 2 are stretching factors (we 
have chosen si = 0.98 and s 2 = 6.5). We then per¬ 
form the same, Re = 10, Poiseuillc flow simulation with 
the three above different grid arrangments with the FV 
method, and for completeness we also include the results 
obtained on a uniform grid by both the FV and the ST 
method. In order to have a better understanding on the 
accuracy of the methods this time we change the num¬ 
ber of grid points (N) while keeping fixed the physical 
domain size L = 64. In other words what we vary here 
is the average grid spacing (Ax) = L/N. 

Figure [5] reports the results of the described test. Three 
sources of inaccuracy lead to the overall error behaviour 
observed here. The asymptotic behaviours are respec¬ 
tively dominated by the spatial accuracy error E/\ x and 
by the finite Mach correction Emu , while the transition 
between these two regimes is also affected by the relax¬ 
ation time error E T . At small resolutions (small values 
of N or equivalently, values of Ax > 1 in our numerical 
experiment) the spatial discretisation error of ST method 
goes as E/\ x ~ 0(Ax 2 ) while the one of the uniform-grid 
FV method behave roughly as E/\ x ~ 0(Ax 3 ) (due to 
the QUICK flux computation). In the opposite limit 
(large V, or equivalently Ax < 1) the compressibility 
corrections comes into play. This effect which goes 
as Emo, ~ 0(Ma 2 ) has a constant behaviour, both 
in the ST and FV method, because in this numerical 
experiment U is kept fixed. This explain the observed 
plateau in the same figure. At the crossover between the 
two regimes, around the value Ax ~ 1 it also happens 
that At ~ 2 t and this corresponds to the range of best 
accuracy on the viscosity term (minimal E T ) (same as 
in figure [3]). In conclusion, there exist an optimum value 
of N linked to the relaxation time error E r for which 
the error is minimum, this happens both for the FV 
and ST algorithms, both on uniform and on stretched 
grids. The ST method (which can only be based on 
a uniform grid) performs better than the uniform-grid 
FV implementation for almost all the values of A, fur¬ 
thermore its absolute accuracy is the highest. However, 
the situation becomes interesting when the non-uniform 
grid FV method is employed. There we notice that 
one can get the same accuracy of the ST algorithm but 
with a smaller amount of grid points. For instance, in 
the case of the hyberbolic-tangent grid with N = 11 
one can obtain the same accuracy as the ST algorithm 


with N ~ 46. This leads to a saving in memory and 
potentially in computational costs. In conclusion, the 
reduction in memory occupation at comparable accuracy 
seems to be the main benefit on can get from employing 
the wall-refined FV method rather than the standard ST. 



FIG. 5. Relative error on the Re = 10 velocity Poiseuille 
flow profile at changing the number of grid points ( N ) and 
keeping fixed the domain size L = 64. Data are traced for 
FV with uniform grid, with Chebychev points and hyperbolic 
tangent grid refinement, and finally for the ST algorithm with 
uniform grid. 

However, a situation that often occur in the simula¬ 
tion practice is that one wants to use all the available 
memory of a computer and using an algorithm with the 
best possible accuracy. The interesting question is then: 
How can we increase the accuracy at comparable mem¬ 
ory costs? Let’s imagine one wants to perform again 
the same Poiseuille simulation at Re = 10 but wants to 
reach a higher level of accuracy (with accuracy defined 
in the sense of L 2 norm). One new possibility is to ad¬ 
just L and U max in a way that the averaged grid spacing 
(Ax) = L/N, with N left unchanged, is the one that of¬ 
fers the best accuracy performance for a given grid. For 
the above case of the hyberbolic-tangent grid this would 
be around (Ax) = L/N = 64/11 = 5.8. The result of 
this novel Poiseuille flow test is shown in figure [6] We 
can see that when N = 64 the best choice is to adopt 
a grid with tanh or Chebychev spacing and with (Ax) 
much larger than unit. The optimum is here reached 

















9 


when (Ax) — 20, this produce an increase in accuracy of 
a factor greater than 100 compared to the case of a simu¬ 
lation with the ST algorithm (and this independently of 
the value of (Ax) chosen for the ST method). 



FIG. 6. Relative error on the Re = 10 velocity Poiseuille flow 
profile at changing simultaneously the domain size L and the 
forcing amplitude, but keeping constant the number of grid 
points N = 64. Here the data are plotted versus the aver¬ 
age grid spacing (Ax) = L/N. Data are traced for FV with 
uniform grid, with Chebychev points, hyperbolic tangent and 
hyperbolic sine grid refinement. The corresponding relative 
error for the ST algorithm with uniform grid is also traced. 


IV. PERFORMANCE EVALUATION 



FIG. 7. Maximum allowed time step in the decaying Kol¬ 
mogorov flow by using eq. © (FV Euler), eq. Q (FV 
Semi-implicit) and eq.© (FV Semi-implicit + Heun). The 
very same result is obtained in the steady Poiseuille flow at 
Re = 10. The horizontal dashed line represents the stan¬ 
dard choice of the time-step for the ST implementation, i.e., 
At = 1 independently of r. 


seen for the simple Poiseuille flow this bring further sav¬ 
ing in terms of computational costs as compared to the 
ST method. 

Finally we shall mention that memory occupation is also 
part of the performance of an algorithm: According to 
our estimate FV in the present formulation needs twice 
more memory allocation as compared to the ST. 


From a computational point of view the FV algorithm 
has more operations per time step than the ST algorithm. 
This comes from the fact that while the streaming pro¬ 
cess can be implemented simply as a shift in memory 
the computation of the flux term in FV involves many 
arithmetics operations. According to our measurement 
the present FV algorithm is about 8-10 times computa¬ 
tionally more expensive than ST algorithm per time step. 
However, as discussed in section III B 21 differently from 
the ST algorithm, in the FV the time-step At is a func¬ 
tion of t. The functional relation linking the maximum 
time-step to r for the proposed time-discretizations can 
be measured and it is reported in figure [7] We observe 
that the method based on ©, semi-implicit integration 
in time of the collision term plus a trapezoidal correc¬ 
tion for the advection is superior to the others. In par¬ 
ticular, for this method A t max > 1 for r > 1/8, that 
is to say that the time-step can be larger than the one 
used in the ST method (which is bounded to the value 
1 for Ax = 1). The most advantageous case occurs for 
r ~ 1/2 - which as we have shown above is also the best 
condition for accuracy - in that case A t max ~ 1.7. This 
reduces the ratio of the computational cost FV/ST to a 
factor 5-6. We note that all this reasoning did not take 
into account the effect of non-uniform grids. As we have 


V. BENCHMARK IN A COMPLEX FLOW: 

HIGH RAYLEIGH NUMBER THERMAL 
CONVECTION 

Several LB Finite-Volume methods proposed in the past 
have been tested just on laminar flows as proof of princi¬ 
ple of the proposed algorithms. Few exceptions exist in 
the literature in which the FV method have been bench- 
marked on much more complex, three-dimensional, de¬ 
veloped turbulent flows. One of such exception is the 
model proposed by Amati et al. [37}, which was probed 
in a three-dimensional plane turbulent channel flow. In 
such case however, the grid wall refinement was based on 
a very simple structure of halved-grid spacing near the 
walls and the accuracy of the method turned out to be not 
satisfactory (the computed mean-velocity profile could 
not properly reproduce the log-law of the wall). 

In this section the proposed Lattice Boltzmann FV algo¬ 
rithm is tested to simulate a complex three-dimensional 
statistically steady turbulent flow. Our choice is here for 
the well-studied flow in the Rayleigh-Benard (RB) cell, 
the prototype of thermal convection driven system [50 ). 
The RB set-up considered in this study deals with a cu¬ 
bic domain (of height H and equal lateral sizes L); it has 
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periodic BC on the lateral walls, while on the horizontal 
walls no-slip and isothermal conditions are imposed. In 
this system, the fluid is heated from below and as such 
(when the heating is large enough and a small pertur¬ 
bation is introduced in the system) an instability arises 
and brings the system into convective condition. The 
dimensionless control parameters are the Rayleigh ( Ra ) 
and Prandtl (Pr) numbers and aspect-ratio Ar = L/H 

M- 

For the LB simulation we use a double population ap¬ 
proach T8|. This means that beside equation m we 
integrate an analogous equation for the distributions 
9 a' 

+ c Q ■ Vg a = — (gl q - g a ) a = 0,...,N pop (19) 

dt Tg 

with equilibrium function g® 9 = (T/p) /® 9 where the 
macroscopic temperature is computed as T = T, a g<x and 
the thermal diffusivity corresponds to n = r g c 2 s . Fur¬ 
thermore, in the equation for f a the forcing term F a is 
assigned in order to model the buoyancy force as rep¬ 
resented in the Boussinesq approximation. In physi¬ 
cal space the added buoyant acceleration has the form 
a = —j3(T — To)g where /3 is the volume thermal expan¬ 
sion coefficient and To is a reference temperature taken 
here as the mean temperature between the ones at the 
top and bottom plates. 

In order to validate the double population approach also 
for the FV method, we first address a rather elemen¬ 
tary simulation in steady convective laminar conditions, 
adapting it from a test case already conducted for the ST 
algorithm in Ref. [3- The system is two-dimensional 
(2D) with control parameters fixed at Ra = 10 4 , Pr = 1 
and Ar = 2.02. The fluid is initially at rest (u = 0), while 
the temperature field is initialised by a linear conductive 
profile, T c (z) = — AT(z/H + 1/2), plus a small perturba¬ 
tion (of order O(10~ 2 )AT) breaking the left right symme¬ 
try. Given the weak, but not negligible, compressibility 
of the simulated flow the initial density stratification due 
to gravity should be also taken into account. This avoids 
the generation of pressure waves at the startup of the 
simulation. We do it via the barometric equation, this 
leads to p{z) = poexp (— cj 2 f3g fj T c (z')dz'), where p o is 
a reference density value taken at temperature To. Note 
that in a 2D system, in order not to suppress the linear 
hydrodynamic instability, the cell aspect ratio (Ar) must 
be slightly larger than 2ir/k c (where k c = 3.117 is the 
wave vector of the most unstable linear mode) jl£j]. In¬ 
deed, when Ar = 2.02 the initial perturbation produce 
an immediate kinetic energy growth and a steady convec¬ 
tive flow pattern establishes. The dimensionless heat flux 
(or Nusselt number Nu) goes from the conductive unit 
value up to around Nu ~ 2.66, see & In figure [8] we 
report the results of simulations conducted with the two 
LB algorithms. We can observe (Fig. [5^) that the tem¬ 
poral dynamics of the dimensionless global heat flux, i.e. 
the Nusselt number Nu(t ), is identical for the two simu¬ 
lations, furthermore they both agree with the analytical 



(a) 



x/H 

(b) 

FIG. 8. Comparison of streaming (ST) and finite-volume(FV) 
LB algorithms in a simulation of the Rayleigh-Benard system 
in steady convective state. The system is two-dimensional, 
and characterised by the control parameters value Ra = 10 4 , 
Pr = 1 and Ar = 2.02. (a) Temporal dynamics of dimen¬ 
sionless global heat flux (Nusselt number Nu(t)) as a func¬ 
tion of time, in dissipative time units to = H 2 /k. The Nu 
steady state value is compared to a value linearly interpo¬ 
lated from Clever and Busse calculations In the in¬ 

set, the grid convergence study displaying the absolute er¬ 
ror of the measured Nusselt number (Nu) with respect to 
Clever and Busse (Nuth) vs. number of grid points in the 
y-direction of a 2D grid, (b) Comparison of temperature iso¬ 
lines in the asymptotic steady state. Levels are taken at values 
T n = To ± n AT /8, with n = 0,1,2, 3. 

asymptotic value given by Clever and Busse (HU . A grid 
convergence study, performed by increasing the number 
of grid points of the same factor in each cartesian direc¬ 
tion, shows that the absolute error on Nusselt number 
as compared to the Clever and Busse value reaches the 
second decimal digit (inset of Fig. [5^). The isocontours 
lines for the temperature field (Fig. [8}j) further display 
the excellent agreement between the two LB algorithms. 
The test exhibits not only the good quality of the present 
FV method but also its consistency with the standard ST 
method also for transient (i.e. time-dependent) dynam¬ 
ics. 
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periodic 



FIG. 9. Cartoon of the three-dimensional Rayleigh-Benard 
system 


z/H 



FIG. 10. The mean temperature profile Tm(z) - averaged 
over time and horizontal planes - as a function of the height 
2 in the cell. To better appreciate the agreement between 
different simulation methods we show here a close-up view of 
the profiles in lower/upper 10% of the cell. 


We then move forward to a more complex case. In par¬ 
ticular, we compare our results with the ones obtained 
by Kunnen et al. [52| for a three-dimensional (3D) 
simulation of a RB system (Fig. [9J characterised by: 
Ra = 2.5 • 10 6 , Pr = 1 and Ar = 2 (see also fEU). In 
this condition the 3D system dynamics is already highly 
chaotic (or moderately turbulent). In [52;| the authors 
employed a direct numerical simulation based on a stag¬ 
gered finite-difference discretisation of the Navier-Stokes 
- Boussinesq equation system. The grid they adopted 
has size ( N x , N y , N z ) = (128,128, 64), it is uniform in the 
horizontal directions and has a sink -type refinement (the 
same as in (1161) ) in the vertical direction. Our benchmark 
is as follow, we perform two series of simulations, one 
with the ST method and the other with the FV approach, 
the dimensionless parameters for the two cases are the 
same as the ones of Kunnen et al ., as well as the number 
of grid points per direction. However, while the ST uses 


a uniform grid in the FV case we use exactly the same 



FIG. 11. Root-mean-square temperature profiles T rms , aver¬ 
aged over time and horizontal planes, as a function of the cell 
height 2 up to the cell center height. In the inset, a zoomed 
in view of the lower 10% of the cell. 


H 

L 

At 

T 

Tg 

AT 

P 


g 

ttot 

FV 640 

1280 

4 

0.5 

0.5 

2 

1.325 • 

10 -4 

1 

1.28 • 10 6 

ST 64 

128 

1 

0.05 

0.05 

2 

1.325 • 

10 -3 

1 

2.48 • 10 s 


TABLE I. Parameter values used for the RB simulations at 
Ra = 2.5 ■ 10 6 , Pr = 1 and Ar = 2: height ( H ) and width (L) 
of the cell, time step amplitude At, fluid (r) and temperature 
(t 9 ) relaxation times, temperature gap across the cell AT, 
thermal expansion coefficient value /? and intensity of gravity 
g. ttot is the total simulation time in numerical units. 


grid as the one adopted in the finite-difference simulation 
1521 . The table Q] reports the numerical values of the pa¬ 
rameters adopted for the two LB simulations. Note that 
the large scale velocity U which is roughly proportional 
to the so called free-fall velocity, i.e. U ~ \J(3gATH is 
the same in both simulations. It is a good practice in 
LB simulations to always keep control of the large-scale 
velocity in order to prevent it to take too large values: it 
is worth reminding that in order to reproduce the incom¬ 
pressible fluid-dynamics the condition U < C 1 is required 
(a commonly accepted rule of thumb in LB practice is 
U ~ 0.1). In order to reach a good convergence of the 
statistical observables in the system the RB simulations 
are carried on for a total time (ttot) which spans over 
several large eddy turnover times (T). We estimate that 
t to t — 12 T for both FV and ST simulations, with T 
computed from the zero-crossing time value of the auto¬ 
correlation function of the total kinetic energy. 

In the figures [TU] and [TT] we show a comparison of the ver¬ 
tical mean temperature profile ( T m ) (averaged over hor¬ 
izontal planes and time) and of the vertical root-mean- 
square temperature ( T rms ) profile, which are defined as 
follows: 
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FIG. 12. Mean turbulent kinetic energy k m (z), averaged over 
time and horizontal planes, as a function of the cell height z, 
and close-up view around the BL peak value (inset). 

We find good agreement among all the three types of 
simulations. Furthermore, we observe that when the 
thickness of the boundary layer At is defined by the 
so called slope definition, At = AT(2 d z T m (z) | z==0 ) _1 
(see figure m both the Kunnen et al. data and the 
FV ones have about 10 points in the thermal boundary 
layer (BL), while the ST despite its remarkable agree¬ 
ment with the other methods, has only 3 points in the 
BL. Small systematic differences can be seen on the 
vertical profile of the mean (turbulent) kinetic energy, 
k m (z) = tt 0 ] L~ 2 f* tot f 0 L f Q L \u 2 dx dy , reported in 
figure IT2l k m (z) has a slower rate of convergence than 
the temperature variance, this is the reason why small 
residual statistical discrepancies remain present here de¬ 
spite of the large number of turnover times of the sim¬ 
ulation. In order to appreciate more sensible differences 
between the FV and ST simulation one has to address 
either observables involving temperature-velocity corre¬ 
lation or small-scale quantities, which are more sensi¬ 
tive to the spatial resolution of the mesh, particularly 
close to the walls. For this reason in figure [13] we com¬ 
pare the time averaged quota-dependent Nusselt number 
Nu(z) = nd z (T) + (u z T)/(kA/H) (where for short (...) 
denotes time and space horizontal averages) and the so 
called Bolgiano length Lb(z) = (/? 3) -3 ^ 2 (e u ) 5 / 4 (eT) -3 ^ 4 
(where e u and £t are respectively the velocity and tem¬ 
perature dissipation rates) (54|. As far as Nusselt num¬ 
ber is concerned, despite a very close mean value, we see 
important differences at the wall. This is due to the com¬ 
bined effect of the boundary conditions and the gradient 


computations in post-processing the data. The tempera¬ 
ture gradient computations involved a second order cen¬ 
tral finite difference scheme at the bulk nodes and a sec¬ 
ond order forward/backward finite difference scheme at 
the boundary nodes. These schemes have been applied 
to both FV and ST algorithms. The FV method exhibit 
wall oscillations which are a factor 10 smaller than the 
ones seen for the ST method, making more reliable the 
total heat flux estimate. Furthermore, in the Lb{z) mea¬ 
sure we observe near a 50% discrepancy at the wall and 
a smaller but non negligible difference in the bulk of the 
cell. Clearly a wall-clustered grid is needed to resolve 
observables built on sharp temperature and velocity gra¬ 
dients. 

We now would like to address the matter of determining 
which LB method is computationally more convenient. 
The choice to have the same U in the FV and ST sim¬ 
ulations has an implication on the determination of the 
large-eddy turnover time and therefore the total number 
of time steps needed to perform a simulation of equiva¬ 
lent physical time-span. The reasoning is as follow: the 
large turnover time goes as T ~ H/U therefore on the to¬ 
tal number of time-steps M for a simulation that should 
span a time T scale as M ~ H/At. It follows that the FV 
simulation will need in this case a number of time steps 
larger by a factor 10/4 as compared to the ST one (see 
table |T]). Since the FV is more expensive than ST by a 
factor 8—10 per time step, we get that the added compu¬ 
tational cost of the FV method is up to ~ 25 larger than 
the ST method. However, such an increase in compu¬ 
tational cost shall be properly weighted by the enhance¬ 
ment in the spatial resolution due to the wall stretched 
grid. An univocal guideline is not available in this con¬ 
text. A commonly employed criterion in the numerics of 
bounded flows is to count the number of grid nodes in 
the BL (another, although less restrictive, rule would be 
to take into account the distance of the first collocation 
point from the wall). Here, if we adopt such a criterion 
the ratio is in favour of the FV method over ST by a 
factor 10/3, that means that we need approximately 3 - 
4 more nodes in the ST simulation. However if we want 
also keep the same aspect-ratio of the simulation domain, 
and since ST is bounded to cubic grids, such an increase 
of the resolution shall be applied to every cartesian di¬ 
rection, which makes the FV grid advantage greater of 
a factor of (10/3) 3 ~ 37. In conclusion, in a ST RB 
simulation we need 37 times more computation nodes to 
perform a simulation with comparable resolution of the 
FV method. By combining the above estimates, we see 
that a simulation of same physical time-span and same 
boundary layer resolution, is about 1 — 25/37 ~ 33% less 
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FIG. 13. (a) The time average Nusselt number Nu(z) as a 
function of the height in the cell, up to half cell. Note that in 
a ideal RB system this quantity should be constant, however 
in numerics, often due to the effect of BC implementations, 
small fluctuations are observed through the cell. It is here 
evident the higher quality of the FV method, (b) Bolgiano 
length, Lb(z), averaged over time and horizontal planes, as a 
function of the cell height z. Note the discrepancies both at 
the wall and in the cell bulk. 


expensive for the FV method than the ST one. In sum¬ 
mary, even if the cost per unit physical time in a FV 
simulation is higher than ST, when a criterion for the 
minimal spatial resolution (particularly near walls or ob¬ 
stacles and in a three dimensional geometry) is chosen, 
the FV method becomes advantageous. 


Finally, we have performed RB simulations at increas¬ 
ingly higher Ra numbers (Ra = 10 s , 10 9 ). All this sim¬ 
ulations have around 10 grid points in the thermal BL 
as shown in the figure [TU No numerical instabilities 
were noticeable as Ra was increased, demonstrating that 
the FV algorithm can deal with turbulence at highly 
Rayleigh number conditions. 


FIG. 14. Mean vertical temperature profile, Tm(z), close to 
the upper and lower plates in the Rayleigh-Benard system 
at Ra = 2.5 • 10 6 (o) , Ra = 10 8 (x) and Ra = 10® (□). 
The Prandl number is Pr = 1 and aspect ratio Ar = 2. 
The thickness of slope boundary layers is indicated by the 
dotted lines. Results were obtained by the FV algorithm at 
resolution 64 x 128 2 , 128 3 and 256 3 respectively. 


VI. CONCLUSIONS 


In this paper we have presented a new finite volume al¬ 
gorithm for the Lattice Boltzmann equation. The new 
method has been validated through a systematic compar¬ 
ison with the standard streaming LB approach, by means 
of test case simulations in laminar as well as in unsteady 
and turbulent flows with heat transfer. The tests have 
shown that the FV has the same order of incompress¬ 
ibility accuracy as the ST algorithm and an improved 
spatial accuracy (third compared to second order), it has 
however a much more elevated computational costs that 
we estimate to be around 8-10 times per time-step. One 
notable advantage of the FV method is the possibility 
to adopt stretched rectilinear grids, which makes it suit¬ 
able for the simulation of turbulent bounded flows. Tak¬ 
ing this into account, i.e. taking into consideration (for 
instance) the minimal number of collocation points re¬ 
quired in a boundary layer for a proper simulation, the 
FV algorithm surpasses the ST method. 

A number of improvement can still be made on the pro¬ 
posed algorithm. First, the boundary conditions can be 
improved. We have noticed that in strongly sheared 
flows, such as channel flow turbulence, some spurious 
oscillations at the boundaries can destabilise the simu¬ 
lations. Second, in stretched cartesian grids the num¬ 
ber of interpolant coefficients to be stored is considerably 
more limited than for other cases of structured grids. In 
the former case the treatment of advection can be fur¬ 
ther optimised compared to the one used in the present 
study, allowing for some extra saving in computational 
costs. 

The FV discretisation proposed in this work builds on 





















14 


the standard streaming algorithm, as a consequence the 
two algorithms do not differ much. The major differ¬ 
ence is of course the way in which the advection is com¬ 
puted. However, their strong resemblance may be useful 
in the development of simulation codes that combine the 
two methods. Therefore, one possible development of the 
present study is to set-up simulations that use the more 
efficient ST method in flow domain regions where a fine 
grid is needed, and the FV method in regions where a 
more coarse grid will suffice (among others, a typical ap¬ 
plication could be the simulation of atmospheric bound¬ 
ary layer with ST at ground level and FV on the upper 
residual layer). It is a prospect that we plan to explore 
in a future work. 
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APPENDIX A. MIXED SEMI-IMPLICIT FINITE 
VOLUME LB SCHEME 


In the following we briefly show how to derive the dis- 
cretised form 0 from the finite volume LB equation 0. 
First the discretisation in time is applied by adopting 
a mixed approach: while for the collision and forcing 
terms a semi-implicit method is used, for the advection 
a simple explicit Euler is implemented. This leads to the 
form: 


jW) _|_ pW _|_ —(/® 9 ( t+At ) — /^ t+At )) + F^ +At ^ 


( 22 ) 


Note that by bringing on the left-hand-side all the terms 
to be evaluated at (t + At): 


fi t+At) ~ ^( f e a q (t+At) - fL t+At) + tF C‘+ At >) = /W - At^c a ■ rij f/W 1 + ^ (f e a q M - /« + tF® ) (23) 
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One now introduces the redefined distribution function At this point the new relaxation time is introduced f = 


fa = fa — §f{fa q — fa + "tF q ). It shall be noted that, r + 4^, and the equation can be written as 


At 


at equilibrium condition, from the above definition we 
get fa q = fa q ~ and therefore the relation can be 

easily inverted leading to 

At , - - , At 


fa fa "F 


2, + A + (24) 
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fa + Y --(fa q - fa) + — F a 


+ fff ~ fa) + A tF a 


(25) 


The above equation coincide with 0, once the correc¬ 
tion factor (l — Al) is applied to the force intensity F a 
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